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Multi-scale turbulence modeling and maximum 
information principle. Part 4 

L. Tao* 


We explore incompressible homogeneous isotropic turbulence within the (fourth-order 
model) formulation of optimal control and optimization, in contrast to the classical works 
of Proudman and Reid (1954) and Tatsumi (1957), with the intention to fix specially their 
defect of negative energy spectrum values being developed and to examine generally the 
conventional closure schemes. The isotropic forms for the general and spatially degener¬ 
ated fourth order correlations of fluctuating velocity are obtained and the corresponding 
primary dynamical equations are derived. The degenerated fourth order correlation con¬ 
tains four scalar functions Di, i = 1, 2, 3,4, whose determination is the focus of closure. We 
discuss the constraints of equality for these functions as required by the self-consistency of 
the definition of the degenerated. Furthermore, we develop the constraints of inequality 
for the scalar functions based on the application of the Cauchy-Schwarz inequality, the 
non-negativity of the variance of products, and the non-negativity of the turbulent energy 
spectrum. We intend to indicate the difficulty for a conventional scheme to satisfy all 
these constraints. As an alternative, we employ the turbulent energy per unit volume as 
the objective function to be maximized, under the constraints and the dynamical equa¬ 
tions, with the four scalar functions as the control variables, which is a second-order cone 
programming problem. We then treat the asymptotic state solutions at large time and 
focus especially on the sub-model where the third order correlation is taken as the control 
variable, considering the computing resources available. 


1 Introduction 

In this part, we investigate incompressible homogeneous isotropic turbulence within the 
fourth-order model, in contrast to the classical works of Proudman and Reid nm and Tatsumi 
H; Our purpose is to resolve their flaw of negative energy spectrum values being developed ([8], 
i)- We intend to demonstrate that, though the simplest three-dimensional turbulent motion, 
homogeneous isotropic turbulence provides us valuable information on why the conventional 
turbulence modeling schemes have their defect and how statistical modeling should be framed. 
Here, by the conventional, we mean that a closure is on the basis of certain equality relations 
assumed among the correlations involved, especially between the highest order and the lower, 
like the quasi-normal adopted in [10] and [IT] . 
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Since we explore a different closure strategy, we need to develop the isotropic tensor repre¬ 
sentation for the general fourth order correlation of fluctuating velocity Wi{'x.)wj{y)wk{z)wi{z'), 
under the supposed incompressibility, homogeneity, isotropy and the intrinsic symmetries from 
the correlation’s definition. Furthermore, we derive the isotropic tensor representation for the 
spatially degenerated fourth order correlation Wi{'x)wj{x.)wk{y)wi{z), considering that the evo¬ 
lution equation governing the third order correlation Wi{x.)wj{'y)wk{z) involves only such a 
correlation and the computational complexity of the degenerated is much less than that of the 
general. The interrelationship between the two is exploited to reduce the number of the scalar 
functions in the representation for the degenerated, which are denoted as Di, i = 1,2, 3,4. 
As a consequence of this representation, we need to re-derive the dynamical equations of evo¬ 
lution for the two scalar functions contained in the well-known isotropic representation for 
Wi{x)wj{y)wk{z) ([lO], [H]). 

We discuss in detail the issues regarding the constraints. To guarantee the self-consistency 
of the definition of the degenerated correlation, the above-mentioned representation needs to 
satisfy additional constraints under more degenerated conditions of spatial positions, such as 

Wi{y.)wj{yi)wk{y)wi{y) = Wk{y)wi{y)wi{y.)wj{yi) 

and 

Wi{x.)wj{'x.)wk{y^)wi{y) invariant under the permutation of {i,j, k} 

These constraints of equality are expected to be satished through adequate structures of Di. 

There exist several sources of inequality constraints for the correlations. One results from 
the application of the Cauchy-Schwarz inequality to the correlations and the structure functions 
in the physical space. The second is from the requirement of the non-negativity of the variance 
of products. The third comes from the non-negativity of the turbulent energy spectrum. These 
inequalities are, to a large extent, neglected by or even unenforcible within the conventional 
schemes, except the limited implementation of realizability. The satisfactions of these inequal¬ 
ities expectedly impose more restrictions on the structures of in addition to those from 
the constraints of equality. We should mention that the Cauchy-Schwarz inequality and the 
non-negativity of the variance of products are natural parts of the closure strategy, because the 
issue of closure in turbulence arises from the average treatment of the Navier-Stokes equations 
and these inequalities are closely related to the mathematical ensemble average operation. 

The above-mentioned constraints are intrinsic to homogeneous isotropic turbulence, and 
their enforcement poses a great challenge to a conventional scheme, because such a scheme 
introduces a set of equality relationships to represent the highest order correlations in terms 
of the correlations of lower orders, and these introduced or added may be incompatible with 
the intrinsic ones, as demonstrated by the specific results of [8] and [9]. One strategy to ac¬ 
commodate all these intrinsic constraints is to adopt one objective function to be optimized, 
constrained by the intrinsic constraints and the dynamical equations of evolution for the cor¬ 
relations, with Di as the control variables. Therefore, the turbulence modeling problem is 
converted into an optimal control problem. For homogeneous isotropic turbulence, we tenta¬ 
tively take the turbulent energy per unit volume as the objective function to be maximized, and 
it is shown that this mathematical formulation is a second-order cone programming (SOCP) 
problem when discretized. 

We formulate the problem in both the physical and the Fourier wave-number spaces, the 
latter makes it easy to represent the constraints and the objective function explicitly in terms of 
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the control variables, which offers us the advantage to use the software packages available like 
‘CVX/MOSEK’ ([T], [3]) to hnd numerical solutions. The employment of Fourier transforms 
introduces higher-dimensional integrals, as a negative consequence. 

We notice the challenge faced by the present optimal control strategy, such as the large 
size of numerical simulation resulting from the numerous constraints and finer discretization 
meshes, which places much more demands on computing resources and algorithms. Also, the 
issue of selection and redundancy of the constraints and the issue of uniqueness of the solutions 
are not yet to be addressed. More information will be gathered from the numerical simulation 
under consideration. 

The present report is organized as follows. We construct the mathematical structure for 
homogeneous isotropic turbulence in Section[2l As done conventionally, we start from the 
Navier-Stokes equations, introduce the correlations up to the fourth order and their symmetries, 
and present the equations of evolution for the correlations. It is then followed by the Fourier 
transforms and the isotropic tensor representations of the correlations up to the fourth order in 
the Fourier wave-number space. The primary dynamical equations for the scalar functions are 
derived, the constraints of equality and inequality are established either in the wave-number 
space or in the physical space, and the issue of objective functions is discussed. Some general 
mathematical properties of the resultant SOCP problem are mentioned. In Section[31 we treat 
the asymptotic state solutions at large time, discussing certain possible characteristics and 
scaling. Due to the restriction of computing resources, we focus on the sub-model where the 
third order correlation is taken as the control variable. 


2 Basic Formulation 


Let us consider homogeneous isotropic turbulence in M^. The fluctuation helds of velocity 
taj(x, t) and (scaled) pressure g(x, t) := p(x, t)/p are supposedly governed by the incompressible 
Navier-Stokes equations. 


_ 0 9wi ^ djwjWk) _ I _ d^jwiWk) 

dxk ’ dt dxk dxi dxkdxk ’ dxkdxk dxkdxi 

we can then construct the following equations for the evolution of the multi-point (tensor) 
correlations up to the fourth order. 


d 

dxk 

d 


Wk{x.)wj{y) = 0, 


d 


d 


^^-wk{x.)wj{y)wi{z) = 0 , —Wi{x.)wj{x.)wk{y)wi{z) = 0 , 


d 


d 


^^Wi{x.)wj{y)wk{z)wi{z') = 0, —M;fc(x)g(y) = 0, -^Wk{x.)wi{y)q{z) = 0 (2.2) 


d 


d 


g^w,(x)wj(y) . 

d 


d 


Wi{yi)wk{yi)wj{y) -h —Wj{y)wk{y)wi{yi) 


d 


- —g(y)M;i(x) + z/ 




\ _ 

+ 7^—^ lM;j(x)M;j(y) 


dxkdxk dykdykj 


(2.3) 
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d 


d 


^^Wi{x.)wj{y)wkiz) + —wi{x.)wi{x.)wj{y)wk{z) 


d 


d 


+ ■^Wj{y)wi{y)wi{x.)wk{z) + —Wk{z)wi{z)wi{x.)wj{y) 


d 


d 


d 


q{x.)wj{y)wk{z) - —q{y)wi{x.)wk{z) - —q{z)wi{x.)wj{y) 


dx^ 

( 




+ 


dy 


dxidxi dyidyi dzidzi 


Wi{x.)wj{y)wk{z) 


(2.4) 




dxkXk 


q{yi)wj{y) = 




dxkdxi 


■wi(x)wt(-x.)wj(y) 


(2.5) 


and 




dxixi 


q{x.)wj{y)wk{z) = 




dxmdxi 


Wl{x.)Wm{^)Wj{y)Wkiz) 




-q(x) q(y) = 




-q{x)wk{y)wi{y) 


dykdyk dykdyi 

Following the conventional treatment of homogeneity ([iQ], [H]), we adopt 


( 2 . 6 ) 


(2.7) 


Uij{r) := Wi{x)wj{y) = Wi{0)wj{r), Uijk{r, s) := Wi{x)wj{y)wk{z) = Wi{0)wj{r)wkis), 


U(ij)ki{r,s) := Wi{x)wj{x)wk{y)wi{z) = Wi{0)wj{0)wk{r)wi{s), 


Uijki{r,s,s') := Wi{x)wj{y)wk{z)wi{z') = Wi{0)wj{r)wkis)wi{s'), 


Q(r) := g(x) q{y) = q{0) g(r), Qj{r) := q{x)wj{y) = q{0)wj{r), 

Qjk{r, s) := q{x)wj{y)wkiz) = q{0)wj{r)wk{s) (2.8) 

Here, r := y—x, s := z—x and s' := z'—x. The dependence of the correlations on t is suppressed 
for the sake of brevity. For the fourth order correlation of velocity fluctuations, we include both 
the general Uijki{r,s,s') and the spatially degenerated 17(jj)fc)(r, s), whose consequences are to 
be explored later. The dehnitions of fl2.8l) result in the symmetry properties of 

Uijir) = Uijk{r, s) = Uikj{s, r) = Ujiki-r, s - r) = Ukij{-s, r - s), 

U(ij)ki{r, s) = U(ji)ki{r, s) = U(ij)ik{s, r), U{ij)ki{r, r) = U{ki)ij{-r, -r), 

U(ij)ki{0, r) = U(^ik)ji{0, r), U(^ij)ki{0, 0) invariant under permutation of {i,j, k, 1}, 

Uijki{r, s, s') = Uijikir, s', s) = Uukjis', s, r) = Uikji{s, r, s') = Ujiki{-r, s - r, s' - r) 

= Ukiji{-s, r - s, s' - s) = Uiijk{-s', r - s', s - s'), U{ij)ki{r, s) = Uijki{0, r, s), 

Q(r) = Q{-r), Qjk{r,s) = Qkj{s,r) (2.9) 
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Considering the zero average velocity field, we also impose the inversion symmetry of 

Uij{T) = t/p(-r), Uijk{r, s) = -Uijk{-r, -s), t/(p)fcz(r, s) = U(^ij)ki{-r, -s), 

Uijkiir, s, s') = Uijki{-r, -s, -s'), Q(r) = Q(-r), Qj{r) = -Qj{-r), 

Qjk{r, s) = Qjk{-r, -s) (2.10) 

With the help of (12.91) and fl2.10p . we insert (12.8p into (12.2p throngh (12.7p to get 
J-C/.,(r) = 0. (^J- +J-)f/,„(r.s) = 0. s) = 0. s) = 0. 

■^Qki{r,s) = 0 (2.11) 

OVk 


dt 


U^,{r) 


A 

dvk 


f^ifcjr'(0; r) + 


A 

dvk 


17jh( 0, -r) 






dvi 


dr~ 


dvkdrk 


u^A^) 


( 2 . 12 ) 


d 


dt 


d d 


dri dsi 


d 


d 


—17i,fc(r, s) - ( ^ + ^ ) U^u)jk{r, s) + —[/(.^^(-r, s - r) + —U(ki)ij{-s, r - s) 


(9s« 


d d d d 

= ^Qjfc(r, s) + ^Qjfc(r, s) - ^Qifc(-r, s - r) - —Q*j(-s, r - s) 

UTi USi UTj USk 


2v 


di^ di^ d^ 

dridri dsidsi dridsi ' ’ * 


(2.13) 




(2.14) 


A AVA AA \--(A d\fd d 

dr I dsi) \dr I dsi) ^ \c?r 


+ ds^ [dri + ^2.15) 


and 


<92 <92 

= -wrwrQfcKr,r) 




dvkdri 


(2.16) 


2.1 Fourier Transforms 

It is convenient to reformulate the above mathematical relations with the help of Fourier 
transforms. The treatment converts the partial spatial derivatives into algebraic operations and 
introduces the turbulent energy spectrum in the Fourier wave-number space. The treatment 
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makes it easy to formulate the constraints and the objective function explicitly in terms of the 
control variables which helps to solve the problem numerically, as to be discussed. Its negative 
side is that higher-dimensional integrals are involved in the formulation as indicated by fl2.17p 
below. 

We employ the Fourier transforms, 

= / dk^ij(k) exp(zk-r), Fijfc(r,s)= / dkdl 1) exp[* (k-r1-s)], 

U(ij)ki{r,s) = / dkdl^(y)fcz(k, 1) exp[?(k-r-F 1-s)], 

7 m 3 xR 3 

Uijki{r,s,s') = / dkdldmf/yfcKk, 1, m) exp[? (k-r1-s-h m-s')] , 

Jr^xR^xR^ 

(5(r) = [ dkQ(k) exp(*k-r), Qj(r) = f dkQjCk) exp(*k-r). 




dkd\Qjk(k,l) exp[^ (k-r-|-1-s)] 


(2.17) 


That the correlations in the physical space are real requires that 

(k) = K,(-k), ^*.,(k, 1) = t7,,fc(-k, -1), ^(F),,(k, 1) = %),z(-k, -1), 

t7*,,(k,l,m) = ^,,fcK-k,-l,-m), Q*(k)=g(-k), g*(k) = ^(-k), 

g;,(k,i) = g,,(-k,-i) 

where the superscript * denotes the complex conjugate operation. 

Substituting (I2.17p into (12.91) and (12.101) and combining with (12.181) . we obtain 


(2.18) 


f/.j(k) = q,(k) = (7„(-k) = e* (k), 

f/i,.t(k, 1) = k) = - 1,1) = Uuj(-k - 1, k) = -f/yt(-k, -1) = -%(k, 1), 

^(ji)fcz(k, 1) ^(jj)zfc(l-) k) k, 1) 1), 


dl 


U T 1) 1) U(^kl)ij{, k 1, 1) 


= 0 , 


d\ 


h^(p')A;Z (Ij k) f/p/jp7(l, k) 


= 0 , 


Uijki{k, 1, m) = (/ijikik, m, 1) = Uukjim, 1, k) = Uikjiil, k, m) = Ujiki{-k - 1 - m, 1, m) 

= Ukiji{-k - 1 - m, k, m) = Uiijk{-k - 1 - m, k, 1) = Uijki{-k, -1, -m) = U*jki{k, 1, m), 

g(k) = g(-k) = g*(k), g,(k) = -g,(-k) = -g*(k), 

4fc(k, 1) = Qk,{\, k) = 4fc(-k, -1) = g*,(k, 1) (2.19) 

which indicate that Uij(k), [/(jj)fc;(k, 1), [/^^^(k, 1, m), g(k) and gjj(k, 1) are real and Uijk{k,\) 
and gi(k) are purely imaginary. Next, substitution of (I2.17p into (12.lip through (I2.16p results 
in 

kk h^fcjr(k) 0, i^kk T Ik) Ukjiik, 1) 0, kj Ukjiik, 1) 0, kk U(ij)kiik, 1) 0, 
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(fcj +/j+ mi) 1, m) = 0, 1, m) = 0, /cfcQfc(k)=0, kkQki(i^,\) = 0 (2.20) 


Q(k) = -^ [ d\Qki{k-\,\) 
kp Jrs 

(2.21) 

4(k) = -w / 

|K| 7r3 

(2.22) 


(2.23) 


d (* 

—Uij{k) + 2i/|k|^t7ij(k) = ikiQjik) -ikjQii-k) +ikk d\ (Uijk{k,l) - ^jifc(-k, 1)) 

at 7r3 V / 

(2.24) 


and 

^^Uijkik, 1) + z/ (|k|2 + |l|2 + |k + 1 | 2 ) .,(k, 1) 

z (^ki Z^) ^j/c(k, 1) % kj Qik(^ k 1,1) "i Ik Qij(^ k 1, k) 

+ i{ki + k) f/(iz)jfc(k, \)-iki Ui^ji)ik{~k - 1,1) - *(-k - 1, k) (2.25) 

In the Fourier wave-number space, the general and the degenerated fourth order correlations 
are related through 

U(ij)kiik,l) = [ dmUijkiim, k,l) (2.26) 

7r3 

following from the link in fl2.9p . whose consequences will be explored in Subsection l2.31 Equa¬ 
tions (12.25p . (I2.23P and (I2.26p indicate that the degenerated U{ij)kl plays the role of an inter¬ 
mediate variable between the rate of change of Utjk and the general Ujjkl, if the latter is taken 
as the control variable. 

2.2 Isotropy 

According to HOI, it is sufficient to formulate the isotropic forms of the tensor correlations 
in the Fourier wave-number space under the supposed isotropy, along with the constraints of 
fl2.19p and fl2.20l) . We have the well-known 

4(k) = 0 (2.27) 

U.A^) = I Ua(k) (2.28) 
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and from [ 10 ], 


^ (m) (k) 

^ 0 ^mn ^p ^2 ~t~ ^np ^m ^^2(^5 ^) 

+ ^pm ^n ^ 2(^5 k, /)^, k + 1 + m = 0 (2.29) 

where Ajm(m) := <5im — rriimm/'rn^ is the second order tensor to enforce the divergence-free 
condition from the incompressibility, the two scalar functions contained in fl2.29l) are real and 
possess the symmetry properties of 


Gi{m,k,l) = —Gi{m,l,k) = —Gi{k,m,l), G 2 {m, k,l) = —G2{1, k,m), k-|-1-|-m = 0 

(2.30) 

For the fourth order correlation within the fourth-order model, we encounter two possible 
representations, the general Uijki(k,l,m) and the degenerated t/(jj)fc/(k, 1), and we need to 
derive their isotropic forms, considering that we explore a closure scheme different from the 
conventional, such as the quasi-normal or its variants. Obviously, equations fl 2 . 2 ip through 
fl2.25p involve the degenerated but not the general, which lends one basis to employ only 
f/(jj)fc;(k, 1). However, the isotropic form of [/^^^(k, 1, m) provides certain restrictions on the 
structure of t/(jj)fci(k, 1 ) and vice versa, and it also offers certain clarities on the issue of closure, 
as to be shown. 

We outline the major steps to determine the isotropic form of Uijki{^, 1, m) below. 

We hnd first a primary general fourth-order tensor function $jjfc/(k, 1, m) of wave-number 
vectors k, 1 and m, which contains the arbitrary scalar functions of the invariants of the wave- 
numbers, /c, /, m, |k-|-l|, |l-(-m| and |m-(-k|. Next, we need to impose the symmetry properties 
of Uijki(k, 1, m) listed in fl2.19p to restrict the structure of the related algebraic operations 
are quite complicated since ^ijki contains a very large number of terms. To avoid this complicity, 
we take a different route. Specihcally, considering that Uijki needs to satisfy the divergence-free 
conditions of fl 2 . 20 l) we employ the projection [101 

Uijkl{K 1, m) = Aji(n) Ajj(k) Axk(l) Au(m) <kijkt(k, 1, m), n = k -f 1 -Km (2.31) 

The multiplication by A/j(n) • ■ ■ AL;(m) eliminates many of the terms and the scalar functions 
contained in <kijki to reduce the above expression to 

GjjKL(k, 1, m) = A/i(n) Ajj(k) Axk(l) Au(m) ^^^(k, 1, m), n = k -h 1 -h m (2.32) 

We then impose the symmetry properties of U^ki directly on the reduced ^[jki- The resultant 
Uijkl of fl2.32p satishes these symmetry properties too, since it can be verihed that fl2.32l) and 
fl2.19p lead to 


A/i(n) ■ ■ ■ Aiz(m) <Fhfc,(k, 1, m) = A/i(n) ■ ■ ■ Ai,(m) <Fh;fc(k, m, 1) 

= A/i(n) • • ■ Au(m) $',^^ ( 111 ,1, k) = A/i(n) ■ • • Au(m) ^'^^((l, k, m) 

= A/i(n) • • ■ Auim) ^^^(-k - 1 - m, 1, m) = A 7 i(n) • ■ ■ ALz(m) ^'kiji{-k - 1 - m, k, m) 

8 


















= A/i(n) ■ ■ ■ Auim) - 1 - m, k, 1) = A/i(n) ■ ■ ■ Auim) $'^fci(-k, -1, -m) (2.33) 

Each implementation of the symmetry properties on results in a fourth-order tensor equa¬ 
tion, a summation of the elements like kiljkkh with their coefficients composed of the scalar 
functions in That is, the equation is a linear algebraic equality of multivariate polynomi¬ 

als of wave-number vector components. For the sake of simplicity, we take the standard form of 
multivariate polynomials as the basis and use its linear independence to infer the values of the 
coefhcients involved. To make the derivation procedure simpler, we start with those symmetries 
in fl2.19p containing the argument — k — 1 — m and update the form of sequentially. 

We may also follow the procedure suggested in cni, by applying the symmetries directly 
to Uijkl of 02.321) so as to determine Here, we face the challenge to deal with the 

large number of terms due to the expansion of A/j(n) • • ■ Aj;,i(m). Whether the two approaches 
produce the same result is yet to be resolved. 

With a lengthy procedure, we obtain hnally 

^/jKL(k, 1, m) = A/i(n) Ajj(k) Axfc(l) Au{m) 

X 6ki D(k, 1, m) -F 6ik Sji D{1, k, m) -F du 6jk D{m, 1, k) j , n = k -F 1 -F m 

(2.34) 

Here, the scalar function Zi)(k, 1, m) possesses the symmetry properties of 

D(k, 1, m) = T)(k, m, 1) = T)(m, k, —k — 1 — m) = D(\, k, —k — 1 — m) (2.35) 


with (k, 1, m) denoting the relevant invariants, {k, I, m, |k -I- 1|, |1 -|- m|, |m -|- k|). These 
invariants will be restricted further to /c, /, m and |k -|- 1 -I- m| by 02.261) . as to be demonstrated 
in Subsection 12.31 

In the degenerated case of U{ij)ki (k,l), we follow the idea similar to the above, enforcing the 
associated constraints of the non-integral form in 02.191) and in 02.201) to infer that 


U{ij)kl{K 1 ) 

= AxkO^) ^iJ ^ki Di{k, /, |k -|-1|) -|- {6jk 5ji + 6ji Sjk) D2{k, I, |k -|-1|) 

-h 6ij Ik ki D^ik, /, |k -h 1|) {h kj + kj Ij) Ik h D4{k, /, |k -h 1|) 

-h ki kj Ik ki D^{k, /, |k 1|) -h h Ij h h D^{1, k, |k -F 1|) (2.36) 

where the scalar functions have the symmetry properties of 

Di{k,I, \k + l\) = Di{l,k, \l + k\), i = 1,2,3,4 (2.37) 

Next, we insert 02.341) and 02.36P into 02.26P to obtain 

Dg = Di (2.38) 

whose derivation will be given in Subsection 12.31 Consequently, 02.36P reduces to 


U{IJ)KL{k, 1) 
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(2.39) 


— ‘^KkiX) ^Lz(l) ^ij ^ki -Di(/c, /, |k + 1|) + {5ik Sji + Sji Sjk) D2{k, I, |k + 1|) 

+ ^ij 4 h D^{k, /, |k + 1|) + {ki + Ij) {kj + Ij) ki D^{k, I, |k + 1|) 
We further constrain Di by imposing the two integral constraints of fl2.19p . 


dm 


dm 


U{ij)KL{m + k, -m) - U(^KL)ijim + k, -m) 


= 0 , 


U(ij)kl{t‘^, k) — U(ik)jl{t‘^, k) 


= 0 


Considering the tensor character of these relations, we analyze them in the special coordinate 
system where k = ( 0 , 0 , k), under k 7 ^ 0, (the case of k = 0 is accounted for through continuity). 
Due to the adoption of this special coordinate system, many integrands involved in the above 
relations are odd functions of mi or m 2 or are invariant under the interchange between mi 
and m 2 , and thus, the associated integrals are trivial. With a lengthy but straight-forward 
procedure to evaluate the above two relations component-wise, we get 


^+00 pm-\-k 

/ dm / (i|m-|-k| 

Jo J\m—k\ 


m 


|m -|- k| 


X 


(|m-|-k|^ — 2fc^ — 3 kmQ) (l — 0^) — 2 |m -|- kp 0^ j Di(|k -|- m|, m, A;) 
-I- ^2 |m -|- k|^ — /c^ — 3 m^ 0^ — 4 kmO^ (l — 0^) k‘^ Dadk -|- m|, m, k) 

— kmQ {kmQ + (l — 0^) D 4 (|k + m\,m,k) = 0 , 


/o 


+00 nm+k |m4-kp —m^ 

dm / dim-I-k| m - j -—- (l — 0^) D 2 (|k -|- ml, m, fc) = 0, 

J\m-k\ |m + k| 


^+00 pm-\-k 

/ dm 
'0 J\m—k\ 


/ d|m -|- k| m |m -I- k| (1 -|- 0^) (Di (m, k, |m -|- k|) — D 2 {m, fc, |m -|- k|) j 
J\m-k\ L ^ ^ 

— kmQ (1 — 0^) D^{m, /c, |m -|- k|) = 0, 

+00 pm-\-k 

dm / d|m-|-k| m |m-|-k| 

J\m—k\ 


X 


0 = 


(1 -I- 0 ^) Di(m, fc, |m -I- k|) — 2 (1 — 0^) D 2 {m, k, |m -|- k|) 

— kmQ (1 — 0^) Dsi^m, k, |m -|- k|) 

— {kmQ + m^) {kmQ + k^) (l — 0^) D 4 (m, fc, |m -|- k|) 
|m -I- kP — m^ — 


= 0 , 


2km 


(2.40) 
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In these integral constraints, we adopt the conventional spherical coordinate system to imple¬ 
ment /jj 3 dm due to its convenience, 

mi = m sind cos 0 , m 2 = m sind sin(/), m 3 = m cos 6 *, 

dm = m"^ sin 6 dm d(j)d6, me(0,-|-oo), 0G[O,27r), d G [0, tt] (2-41) 

followed by the change of variables 0 = cos 6 and 0 —)■ |m -|- k| in order to match with the 
arguments of the scalar functions involved. 


2.3 Relationship Between General and Degenerated Fourth Order 
Correlations 

In this subsection, we give a rather detailed analysis on what invariants should be involved 
in Zi)(k, 1, m), why 02.381) needs to hold, and what possibly additional constraints are present 
for D, by exploiting fully 02.261) . This analysis also helps to simplify the treatment of 02.25P if 
one intends to solve for Uukl and D, since 02.251) involves directly U{ij)kl which acts as an 
intermediate variable. 

Substituting 02.34p and 02.36p into 02.261) . we obtain 


Airfc(k) Ai,i(l) 


5ij 6ki Di{k\ l\ |k + in + 6ji + 6n 6jk) D2{k\ l\ |k + Ip) 


+ ^ij {kk + 4) {h -|- li) D^{k^, |k -|- Ip) 

+ (4 kj + ki Ij) {kk + Ik) {h + k) D^{k\ l\ |k + Ip) 
+ ki kj {kk -\- Ik) {ki k) D^{k^, |k -|- Ip) 


+ lilj{kk + lk){ki + li)D,{l\e,\k + l\^) 

= AKk{k) Ali{1) / dmAji{n)Ajj{m)(6ij6kiD{m,'k,\) + 5ik6jiD{k,m,\) 

Jr3 V 


6ii6jkD{\,k,m] 


n = m 


k + 1 


(2.42) 


Here, we have adopted an equivalent set of (squared) arguments for Di for the sake of conve¬ 
nience in discussion. Equality 02.42p is a tensor equation, and because of this tensor character, 
we analyze the equation in the special coordinate system where 


ki — li — 0 , k -|- 1 — ( 0 , 0 , |k -|- 1 |) 


(2.43) 


This is the case of the concerned k and 1 forming a plane, and the coordinate system is rotated 
such that 02.43P holds in the rotated system; The exceptional case of k and 1 parallel to each 
other or zero can be dealt with through the continuous distributions of the scalar functions 
and a limiting procedure. This special coordinate system helps to simplify the mathematical 
operations and motivates us to adopt the form of {kk -l- 4)(4 + k) in the above equation. One 
consequence of (I2.43p is that both |m -I- k -I- 1| and D{m, k, 1) are even functions of mi, which 
makes many integrals involved in fl2.42p automatically zero. 
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The interrelationship fl2.42p needs to be evaluated component-wise. Firstly, the components 
of IJKL = 2311, 3211 result in 


/ dm A 2 j(n) A 3 j(m) Zl(m, k, 1)-h / dm Ai 2 (n) Ai 3 (m) Zl(k, m, 1)-h dl(l, k, m) 

Jm.3 Jr3 L 

/ dm A 3 j(n) A 2 j(m) Zl(m, k, 1)-h / dm Ai 3 (n) Ai 2 (m) Zl(k, m, 1 )Zl(l, k, m) 

jR3 d]R3 L 

n = m + k + 1 


= 0 , 

= 0 , 
(2.44) 


These two integral equalities supposedly contain only the invariants, k, I and |k-|- 1|, reflecting 
the nature of isotropic turbulence. We have assumed that D depends possibly on the invariants 
of k, 1 and m, 

Zl(k, 1, m) = D[k^, |k -h Ip, |1 -|- mp, |m -f kp) (2-45) 

with 


|1 -|- mp = — 2(^2 m 2 + k^ m 3 — |k -|- 1 | m 3 ), |m -|- kp = k'^ + + 2(^2 m 2 + k^ m 3 ) 

(2.46) 

To guarantee the non-explicit presence of components k 2 and k^ in fl2.44p . we need to combine 
the above two quantities so as to eliminate the components, 

|l-l-mp-|- |m-|-kp = k"^ + P + 2 mP -|- 2 |k -|- 1| m 3 (2.47) 

which in turn can be equivalently replaced by 

|k -I- 1 -I- mp = |k -I- Ip -|- m^ -|- 2 |k -|- 1| m 3 (2.48) 

along with /c^, P and m^. Next, we need to exclude the presence of |k -|- Ip as an independent 
argument in D of 02.451) in order to satisfy the symmetry property 02.351) and to guarantee 
no separate presence of |1 -|- mp and |m -|- kp in the integrands of 02.44p . Otherwise, such 
separate presences would result from the combined effect of both |k-|-lp as one of the invariant 
arguments of 02.451) and the various interchanged positions of k, 1 and m in H in 02.441) . Thus, 
to have the desired invariants, fc, I and |k -|- 1|, present in 02.44p . D depends possibly on the 
invariants as follows, 

Il(m, k, 1 ) = D[\m -|- k -|- Ip, m^, /c^, P) (2.49) 

constrained by 

Zl(|k-|-l-|-mp,fc^,/^, m^) = Zl(|k -|- 1 -|- mp, m^, /^) = D{lP^ |k -|- 1 -|- mp, /^, m^) 

= m^, |k -I- 1 -I-mp, (2.50) 

following from 02.35p . 

Under 02.43P and 02.49p . we have the property that 

|m -|- k - 1 - 1| and D{m, k, 1) are even functions of mi and m 2 and invariant under mi ■H- m 2 

(2.51) 
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This property guarantees the automatic satisfaction of fl2.44l) and is used to simplify the con¬ 
straints of other components below. 

Secondly, the components of JJ = 12,13, 21, 31 and the symmetry properties of D and D 2 , 
(1051) and flQTD . yield 


zl2(fc^/Mk+l^ 



(^Aii(n) A22(m) Ai2(n) Ai2(m) jT)(k, m, 1 ) 
(^Aii(n) A33(m) Ai3(n) Ai3(m)^^(k, m, 1 ) 
(^Aii(n) A33(m) -h Ai3(n) Ai3(m) j D{1, m, k). 


n = m -|- k + 1 
(2.52) 


where the integrations are implemented with the help of fl2.41l) . 

Thirdly, the components of IJKL = 1111,2211,3311 result in 

Di(fc 2 ,/Mk+in 

= / dm Aij(n) Aij(m) Zl(m, k, 1) -|-2 / dm Ai 2 (n) Ai 2 (m) Zl(k, m, 1) 

= / dmAj 3 (n) Aj 3 (m)T)(m, k, 1)-h 2 / dm Ai 3 (n) Ai 3 (m) T)(k, m, 1) 
dR3 dR3 

= / dm Aij(n) Aij(m) dl(m, k, 1 )-|-2 / dm Aii(n) Aii(m) Zl(k, m, 1) 

Jm.3 J]r3 

— 2 dl 2 (/c^,|k -I- 1|^), n = m -I- k + 1 (2.53) 

The last equality indicates that these expressions for Di can also be cast as the expressions for 
D 2 . 

Fourthly, the component of IJKL = 1122 leads to 




Ik + ll 


dm Ai3(n) Ai3(m) - Ai2(n) Ai2(m) m(k, m,l). 


n = m -|- k + 1 


(2.54) 


which provides the integral relationship between D 3 and D. 

Fifthly, we evaluate the equalities under IJKL = 2222, 2322 and employ the linear inde¬ 
pendent properties of { 1 , (^ 2 )^} and { 1 , ^ 3 , (^ 3 )^} to obtain 

Di = D 5 (2.55) 


Finally, the equality from IJKL = 3322 gives 

2 Zl2(A;^/Mk + in + |k + 1|2113(^2,/Mk + 1|2) + |k + 1|4 Zl 4 (fc^/Mk + in 
= 2 [ dm(A 33 (n)A 33 (m) - Ai 3 (n) Ai 3 (m))T)(k,m, 1), n = m-F k-h 1 (2.56) 

JRS ^ ^ 
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which provides the integral relationship between 1)4 and D. 

A few comments are worth to make here. The multiple expressions for D 2 in fl2.52p and 
for Di in fl2.53p represent a set of the integral constraints of equality for the scalar function 
D. Based on the role played by U(ij)kl as the intermediate variable, these constraints may 
also be expectedly derived from the substitution of fl2.26p and fl2.34p into the tensor equation 
of evolution fl2.25p : the component-wise satisfaction and implementation of the latter produces 
these constraints and the dynamical equations (I2.58p and (I2.59p below. These constraints need 
to be imposed explicitly, if Uijkl (or D) is adopted as the control variable in the fourth-order 
model. This adoption poses a great challenge computationally ~ how to handle adequately 
D(k, 1 , m) = Zi)(|k -|- 1 -|- mp, , m?) which results in a larger number of discretized control 
variables and constraints, the integrals of higher dimensions, etc. Additionally, there is a subtle 
issue in the treatment of Uukl (or D) as a control variable, since the evolution equations, 
(I2.2ip through (I2.25p . involve only U(^ij)kl and additional information contained in the solution 
of Uijkl expectedly comes from the related constraints and optimization. Next, the required 
satisfaction of the constraints for D indicates the difficulty to approximate Uijkl appropriately 
in terms of lower order correlations as done in a conventional closure scheme, because of the 
great difference between the information contents in these correlations of various orders and 
various numbers of vectors involved. 


2.4 Primary Dynamical Equations of Evolution 

With the help of the isotropic forms presented in fl2.27p through fl2.29p and fl2.39p and with 
the help of the special coordinate system in which either k = (0, 0, k) or (I2.43p holds, we can 
derive from (I2.24p and (12.251) the following dynamical equations governing the evolution of the 
scalar functions Ukk, Gi and G 2 , 




+ 2uk‘^ ]Ukk{k) 


^+00 rl-\-k 

= 47r / dl d|k + l|/|k + l| (1-02) 

Jo J\l-k\ 


X 


\l-k\ 

G^k + IQ) 

|k + l |2 


, ^ 2Ue^ + kiQ-P ,, ^ 

— ( 2-1-- ) A;G 2 (|k -I-1|, Z, fc) 


(1 - 02 ) A :2 Gi(|k - 1 - 1 |, fc, /) + ( 0 + ] I G2{\k + \\,k, 1) 


k + 1 

, 0 = J 


k-hl |2 - A ;2 


2 kl 


(2.57) 


,A:,0=0 (2.58) 

and 

+ + f + \k + \\^)^G2{l,\k + \\,k) = D2{\k + l\,k,l) - D2{\k + \\,l,k) (2.59) 

These dynamical equations have several interesting features. The first is that equation fl2.57p 
is another version of the known (see Equations (52) and (54) of llQ]); the change of variables. 


(^ + ^(^' + ^' + |k + ir))Gi(|k + 
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0 —|k + 1|, is adopted to be suitable to the arguments of the scalar functions involved. The 
second is that the evolution of Gi has a zero source, the same as Equation (59) of jTD]. The 
third is that the evolution of G 2 is directly affected only by D 2 , the impacts of Di, and 
are through the integral constraints of fl2.40l) and the constraints to be formulated. 

With the help of fl2.37p . we can infer from fl2.59p that 

+ /' + |k + ir))(G2(/, |k + 1|, fc) + G2(|k + l\,k,l) + G2{k, I, |k + 1|)) = 0 (2.60) 

which is equivalent to Equation (62) of [lO]. Both (12.581) and (I2.6np may be useful in the 
analysis of asymptotic state solutions at large time. 


2.5 Constraints 

Within the fourth-order model with Ui^ij)KL, or equivalently Di, i = 1,2, 3,4, as the control 
variables, the following sets of constraints need to be satisfied. 

Firstly, we have the constraints of equality (I2.37p and (I2.40p . which involve only the scalar 
control variables. These constraints are to guarantee the self-consistency of the dehnition of 
the degenerated fourth order correlation introduced in fl2.8p . 

Secondly, to satisfy the constraints of equality fl2.30p . we need to impose them on the initial 
conditions of Gi and G 2 , as required and guaranteed by fl2.58p and (12.591) . The equality of 
fl2.27p is satished automatically, which can be directly verified under k = (0, 0, k). 

Thirdly, the non-negativity of the energy spectrum /c^ Ukk{k)/2 requires that 

Ukk{k) > 0 (2.61) 

This inequality has several consequences similar to those discussed in PART III [I3] . It results in 
the positive semi-definiteness of the second order vorticity correlation cUjCUj (k) defined through 

Uj{^) = ejkiWk,iix.), Ui{x.)uj{y) = / dka;ia;j(k) cos(k ■ r), QiUj(k) = k"^ Uij(k) (2.62) 


It leads to the following inequality for WiUj^k), 


WiUJj{k) 


WiLLlj{k) 


2 


< Uu{k)ujj_Uj_(k), 

% ~ 

2 ^jil ki Ukki^k^ 




dkwiUjik) sin(k ■ r). 


(2.63) 


Furthermore, the corresponding inequalities held automatically in the physical space can be 
established, as done in [13], 


%( 0 ) > 0 , 


^p(r) 






— 2 ^fc/c( 0 ) 


(2.64) 


cOiUJiiO) > 0, (^a;ia;j(r)) < u^u}j_{0), UiUj{r) := -eimn(-jki 


d'^Umki^) 

dvndri 


(2.65) 
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( 2 . 66 ) 


2 dJJ-ihiv^ 

< ?7ij(0) ayJ7(0), M^(r) := 


dri 

which have the characteristic of spatial degeneracy on the right-hand sides of the inequalities. 

Fourthly, besides the above resultant and redundant inequalities, we can also generate the 
constraints of inequality involving only the second order correlation by applying the Cauchy- 
Schwarz inequality, \ah\ < aabb, to structure functions [2] involving only the second order 
correlation such as 

[wi(y) - u!i{x) + a{wi{z') - Wi{z))] [wj{z') - Wj(z) - /3{wj(y) - 'iCj(x))], 

[wi(y) - u!i{x) + a{wi{z') - Wi{z))] [ujj{z') - ujj{z) - - a;j(x))], 

[a;i(y) - a;i(x) a(ui(z') - u:i{z))] [u:j{z') - u:j{z) - (3{u:j{y) - a;j(x))], a > 0, /3 > 0 (2.67) 
which results in the following primary constraints of inequality, 

Uij{s' - r) - Uij{s - r) - Uij{s') + Uij{s) + 2a{Uij{0) - Uij{s' - s)j 

-2^(Uij{0) - 17,,(r)) - a/3(u,,{s' - r) - 17,,(s - r) - t/„(s') + 17,,(s)' '' 

< 4 [f/,i(0) - Uiiir) + (Uii{0) - [/^(s' - s)) 

+ a(Uii{s' - r) - Uii{s - r) - Uii{s') + Uii{s) 

X U,J0) - 7/,,(s' - s) + (t/,,( 0) - f/,,(r)) 


7/„(s' - r) - 7/,,(s - r) - 7/,,(s') + 77„(s 


, *j = ll ,12 


( 2 . 68 ) 


WiUj{s' — r) — WiU)j{s — r) — WiUj{s') -|- WiUj{s) 

+ a/S ( WiUj{s' — r) — WiUj (s — r) — WiUj{s') + WiUj{s 


<4 


Uii{0) - Uii{r) + a‘^(Uii{0) - U^is' - s)j 

+ a(Uij{s' - r) - 74i(s - r) - Uii{s') + Uij{s 
ujjUJj (0) — oJjOJj (s' — s) -f {^j^j (0) “ Ci;,Ci;, (r) j 
-mujujis' - r) - (J]uJj{s - r) - (J]uJ]{s') + uJuJjis 


and 


, *j = ll ,12 


UiUj{s' — r) — UiUj{s — r) — a;ja;,(s') -|- UiUj{s) -\- 2a{uiUj{Qi) — UiUj{s' — s) j 


(2.69) 
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<4 


— 2/9 ^cJiCJj(O) — uJiUJj{r)^ — a(3 (uiUj^s' — r) — UiU)j{s — r) — UiUj{s') + cJiCJj(s) j 
cu^(O) - cu^(r) + ^a)^(O) - uJm{s' - s) 

+ «(cuffs' - r) - (AuJlis - r) - cu^(s') + w^(s' 


X 


ix;ja;j(0) — UjUj{s' — s) + (ujUj{0) — UjUj{r) 


- - r) - - r) - 


, *J = 11,12 


(2.70) 


Here, each inequality involves essentially a scalar function defined in a high-dimensional space 
composed of the components of r, s, s' and their differences and t. Additional inequalities 
can be produced similarly based on other structure functions. We expect that these quadratic 
constraints play a signihcant role to constrain directly the structure of Ukk and indirectly the 
structures of the control variables. 

Fifthly, for the higher order correlations, we apply the Cauchy-Schwarz inequality to the 
correlations of 


Wi{y.)wj{y)wk{ 2 ), Wi{-K)wj{iL)wK{j)wL{'z), g(x)g(y), g(x)M;i(y)M;^(z) (2.71) 


Here, more are to be added like the correlations involving tCfc ;(x), a;j(x), and the structure 
functions involving Wi{y) — tCi(x), a;j(z') — Ui{z), etc. In the derivations of the primary con¬ 
straints of inequality below, we resort to (12.91) and the interchangeability between the space 
vector components {vj, Sj}, j = 1, 2, 3, like {r 2 , S 2 } •H- {ri, Si}. 

1. For Wi{'x)wj{y)wk{z), there is only one independent decomposition. 


Uijkijl'j s) 


< 2 Uii{0) j j(r, r), 77(ii) j j(r, r) > 0 


(2.72) 


which in turn results in the primary constraints as follows. 


U(ii)jjir,r)>0, i = l,j = l,2-, 


Uijk{r,s)\^ = 111, 112,121,123 (2.73) 


2. In the case of Wi{'K)wj{x.)wKiy)wL{z), we need to decompose the quantity properly such 
that the decomposition correlations are resolvable within the fourth-order model. 


U{ij)KL{^^ s) < 0) ^(J£K)ll(^ ~ r, s — r), 

2 

U{ij)KL{^,^) ^ r) f^(ji)LL(s, s) 

whose primary constraints are 
2 

U(ij)KL{^i^) < j(0, 0) — r, s — r). 


(2.74) 
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ijKL = nil, 1112,1122,1123,1211,1212,1213,1233; 

2 

U{ij)KL{'^-i^) < r) f/(j s), 

i]KL = nil, 1112,1122,1123,1211,1212,1213,1221,1223,1233 


3. For g(x)g(y). 


|Q(r)|<Q(0), Q(0)>0 


4. In the case of q{'s)wi{y)wj{7.), we obtain 


Q*j(r, s) < Q(0) - r, s - r), ij = 11, 12 


(2.75) 

(2.76) 

(2.77) 


Sixthly, we formnlate the constraints of ineqnality on the basis of the requirement that the 
variance of products be non-negative, 

(Xr-ZF)^>0, XXYY > (XYf (2.78) 

We apply it to the specihc cases of 

= (^*(x),Wj(y)), (wi(x),a;j(y)), (a;i(x),a;j(y)) 

to generate 

t/(ii)^(r,r) > (^f/jj(r)j , M;i(x)M;i(x)n;j(y)a;j(y) > (^M;ia;j(r)j , 

a;i(x)a;i(x)a;j(y)a;j(y) > (^a;ia;j(r)j , zj = 11,12 (2.79) 

More such inequalities will be produced like setting X = Wk^ii'x), Wi{y) — tCi(x), ujj{z') —u)j{z), 
and so on. 

The constraints of equality and inequality listed above are an integral part of the mathemat¬ 
ical setup to model incompressible homogeneous isotropic turbulence. The issue of how to deal 
with numerous constraints constructed from structure functions and the issue of redundancy 
are not addressed; this subject needs to be explored due to its importance to the geometry 
and size of the domain of feasible solutions and the computational feasibility. Each of the 
above constraints may be rather loose; however, their number is substantial and the very many 
such constraints may form a rather tight restriction on Di. These constraints are either linear 
or quadratically convex, when discretized, as functions of the discretized control variables Di, 
z = 1, 2, 3,4, since fl2.57p through fl2.59p are of linear structures, similar to that argued in [T3] . 

2.6 Expressions for Correlations in Physical Space 

We present here the expressions for the correlations in the physical space which are needed 
for the implementation of the constraints formed in the physical space. 
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To this end, we substitute fl2.2ip . fl2.23p . fl2.28p . fl2.29p and (12.391) into fl2.17p and then 
operate on the resultant expressions with some conventional techniques. The first technique is 
to use 


kj exp(*k-r) 


1 c?exp(zk-r) 
^ dvj 


and the like involving Ij. With the aid of this technique, we can replace all the wave-number 
components of free indexes in the expressions with the corresponding partial derivatives with 
respect to or Sj. For example, we obtain 






dkUkkik) exp(?k-r) -h - 




2 dvidr^ 


dk — Ukk{k) exp{i'k-r) (2.80) 


The second is the adoption of the spherical coordinate system fl2.4ip . k —)■ {k,(j),9) and 1 —)■ 
in order to introduce r = |r| (and/or s = |s|) into the expressions and to simplify 
the integrations analytically. Specifically, under given r and s, we transform/rotate k and 
1 such that k ■ r = r k^ and 1 • s = s /3 (or 1 • k = /c /3 if s is absent and 1 present). The 
third is to integrate 9 (O') and 0(0') analytically, if possible, such as d(p c/0'/(cos(0 — 
0 ')) = 27r Jq''c/ 0 /(cos 0 ) = 47 r JJ'c/ 0 /(cos 0 ) = 47 r <^ 0 /(“ sin 0 ), and employ the change 

of variables 0 = cos 6*, 0' = cos 9', etc. 

Using the techniques and procedure outlined, we can derive from (12.171) the following 


Q(r 

: 87r 


+00 ^+00 ^1 j 2 

2 / / dl I dQ' — 

k^ 


/ dk / 

'0 Jo 


'-1 


x{Di{\k + l\,l,k) r k 


2 , ,2 {k^ + k-\y (k-l)2 (fc2 + k.l)(k-l)(k-l + /2 


Ik-hll 


P 


+ 


Ik-hll 


P 


+ (2 D2i\k + 1|, /, fc) + k^ Zl3(|k + 1|, /, fc) + k^ D^{\k + 1|, /, fc)) 


X k^- 


(/c^ -h k ■ 1)^\ f 2 (k ■ 1)^\ 1 sin(r/c) 


kie', 

|k + l|‘ 

u.,M 

27r r 

r \ 


|k + l| 

, i2 


r- 


p 


rk 


r if- / 1 \ sin(r/c) cos(r/c) . 

5ij J dkk [ sm(r/c)- ^ 2 , 2 - ~i — ) ^kk{k) 

r+oc 

I dk '-' " - 


TiTj 


r‘^k'^ ' rk 

{r‘^k‘^ — 3) sin(r/c) -h 3r/ccos(r/c) 
k 


(2.81) 


(2.82) 


QKL{r,s) 
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X <'Di(fc,/,|k + l|)|k + l|^( -5kl- 
' 2 D 2 




92 k ■ 1 

+ 


k"^ drxdrL P OskSsl k‘^P OtkOsl, 
{k,l,\k + l\) + \k + l\^D3{k,l,\k + \\) + \k + \\^D^{k,l,\k + l\)^ 

^ I cos{rkQ + s/0') (2.83) 


d 


k-1 d 


d 


k-1 d 


dsK k'^ dvK J \drL P dsi 


Uijkir,s) 

p-\-oo p-\-oo pi pi r'^!^ 

= 47r / dk dl de dO' / d^kH'^ 

Jo Jo J-l J-l J-it/2 

d + k • 1 / a d 

dvi |k + ip 9si 


and 


X <( — Gi(m, k, 1) 
+ G2{Tn, I, k) 


d k-1 d 




k-1 d 


dsj kJ dvjj \drk P dsk 


1 1 / d 

dvidrj + ik+iK dri (9sj / \ dvj dsj 


A:2 + k -1 , 

^ d 

pd\d] 

( ^ - 

k-1 d 

|k + l|2F 

^dvi 

dsi J dr,j 

\drk 

P dsk 


+ 02 ( 1 , m, k) 


d + k -1 / a d 

dvi |k + ip \c?rj (9si 

/ 1 (92 1 (92 k -1 (92 

1 jfc + p Qg Qg^ /j2 Qj- drk kJ P dvjdsk 


+ G2{rn, k, 1) 


■ 1 92 1 

Oki + —7^- 7 - 1 “ 




d d 
+ 


P dskdsi |k + 1|2 \drk dsk J \dri dsi 

k-\ + P d f d k-ia 

dsj kJ dvj 


\k + \\'^ P dsk\dri (9s. 


X cos(r/c0 + s/0') 


(2.84) 


U{IJ)KL{r, s) 

^+00 ^+00 pi pi /*'^/2 

=An / dk dl de de' / d^k^P 
Jo Jo J-l J-l J--K/2 


(92 


X <1 Di{P /, |k + 1|) 5ij + ^2 Qr^dvL ' p OskOsl k^ P OtkOsl 


1 

+ 


+ D2 {kp,\k + l\)( 6 iK + ^jr^ 

\ k^ orjoVK 

+ D2{k, /, |k + 1 |) (6jk + 7-7 


-Ds{k,l,\k + l\)Sij 


kJ drjdrx 

d k-1 (9 


5jl + - 

5/l + ^ 


92 k-1 (92 

1 92 


P dsjdsL 
1 (92 

P dsjdsi 

d k-1 (9 


dsK /^2 drx J \drL P dsi 
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d d 
+ 


d k-1 a 


X 


d 


dr I dsj J \drj dsj J \dsK k"^ drx 

k ■ 1 5 II / , ^ 

cos{rkQ + slQ) 


dri P ds L 


Here, 


k-l = kl(e 0 ' - Vl - 02 Vl - 0 '^ 


sin 


k + l|2 = fc2 + /2 + 2k-l 


(2.85) 


( 2 . 86 ) 


are used in 02.831) through 02.851) . 

The integrals can be evaluated numerically with the help of software packages like CUBA 

(i, 0, 0). 


2.7 Objective Function 

Within the supposed incompressibility, homogeneity and isotropy, we have set up a mathemat¬ 
ical structure with the correlations up to the fourth order, without additional approximations 
involved. This structure contains three fundamental elements: One is the primary dynamical 
equations of evolution for Ukk, Gi and G 2 , fl2.57p through fl2.59p . The second is a set of linear 
equality constraints and a set of linear and quadratic inequality constraints listed in Subsection 
fm The third is the scalar functions Di, D 2 , and £> 4 , yet to be determined. 

Due to the linear forms of the dynamical equations fl2.57p through fl2.59p . Ukkit), Gi{t) and 
G 2 {t) can be formally solved in terms of D 2 {t), t < t, under appropriate initial conditions. 
As a consequence, the above-mentioned three elements effectively result in a mathematical 
structure of linear and quadratic constraints intrinsic for Di, D 2 , and D 4 in incompressible 
homogeneous isotropic turbulence. 

It is known from the closure problem of turbulence that this mathematical structure of 
intrinsic constraints allows many feasible solutions for Di to exist, under an initial condition 
given; and a specihc closure tends to reduce the number of solutions. A conventional closure 
scheme, like the quasi-normal, essentially adds another set of equality constraints, and these 
added ones may not be compatible with the intrinsic, as specihcally demonstrated in [S] and 
[9]. Considering the big number and various origins of the intrinsic constraints, we expect 
that any such closures will have the difficulty to satisfy them. This observation leads to the 
option to bypass the conventional: one objective function is introduced which is to be optimized 
under the intrinsic constraints, with Di, D 2 , D 3 and D 4 as the control variables. That is, the 
turbulence modeling problem is transformed into an optimal control problem. 

The question now is what objective function is to be adopted. We select tentatively the 
turbulent energy Ukk{0), based on the following considerations. The hrst is that it reflects the 
accumulative impact of the third and fourth order correlations, as indicated by the dynamical 
equations fl2.57p through fl2.59p . The second is that it is an invariant trace of the second order 
correlation tensor or matrix and it may be used to quantify the spread-out of the probability 
density function of the fluctuations, as argued in PART I mi- The third is that it results in a 
linear (and concave) function of Di, D 2 , and D 4 to be maximized, which have the advantage 
to be more easily dealt with from a computational perspective, unlike any other invariants of 
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the second order correlations. It then follows that 

Ukk{^) = / dk 17fcfc(k) = dvr / dkk^Ukk{k) to be maximized (2.87) 

Jrs Jo 

for the problem of our concern here. The objective is a linear function of the control variables 
Di, D2, D3 and D^, when discretized, with Ukk, Gi and G2 as the primary state variables and 
under the intrinsic constraints listed in Subsection 12.51 The present optimization problem is a 
second-order cone programming problem [7j. 

We mention that there is a thermodynamic setting underlying the Navier-Stokes equations 
fl2.1l) in that ly > 0 for an incompressible Newtonian fluid is imposed on the basis of the second 
law of thermodynamics. Dealing with the statistical average of the hydrodynamic turbulent 
fluctuations, we construct the present formulation from the perspective of information theory 
regarding the objective function; The intrinsic constraints result from the mathematical average 
operation - the supposed isotropy and the applications of the Cauchy-Schwarz inequality and 
the non-negativity of the variance of products to the elementary correlation functions and the 
structure functions; The impact of the Navier-Stokes equations is reflected by the dynamical 
equations. We have not included any constraint built on the basis of specific physical models, 
considering its potential incompatibility with the intrinsic ones. 

2.8 General Properties 

Some mathematical properties of the SOCP problem can be established straight-forwardly. 
We have the known relation of 

(9 ( 9 ^ P ~ 

= 2u—^U,,{v) =-2u d\ieUkk{k) < 0 ( 2 . 88 ) 

at clrkOrk j.=o 

following from fl2.12l) . fl2.17p . fl2.27p . fl2.29p and fl2.6ip . This monotonic decay is expected 
physically due to the lacking of external energy supply and the effect of viscous dissipation. 
Next, we have that 

if {Di, D2, D3, D4, Gi, G2, Ukk} is an optimal solution, 

then A {Di, D 2 , D 3 , D 4 , Gi, G 2 , Ukk}, A G (0,1], is an optimal solution too. (2.89) 

Here, the restriction of the scaling factor A < 1 comes from fl2.79p . 

Thirdly, a known scaling property of the Navier-Stokes equations implies that under the 
transformation of 

|f, k, r, Di, D2, D3, D4, Gi, G2, 

k, r, z/'^Z 92 , z/^Gi, v^G2, i^'^Ukk^ ( 2 . 90 ) 

the kinematic viscosity u will be removed from the resultant SOCP problem whose mathematical 
equations have the same forms as the original under z/ = 1. 

Compared with the formulation of HDl and [H], the present is much more complex mathe¬ 
matically and computationally, resulting from the different closure strategies underlying. 
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3 Asymptotic State Solutions 

To test the SOCP problem above against the experimental and DNS data, we focus on 
the asymptotic states of decaying at large time, which are characterized within the present 
framework by certain conditions to be specihed below. 

First, as pointed out in [10], dynamical equations fl2.58p and fl2.60p can be solved directly 
for Gi(|k + 1|, A;, /, t) and G2{1, |k + 1|, k, t) + G 2 (|k + 1|, A;, /, t) + G 2 {k, I, |k + 1|, A) which decay 
exponentially with time under k 7 ^ 0 or 1 7 ^ 0. Consequently, we take these quantities as 
essentially trivial at large time to characterize the asymptotic states, 

Gi{m,k,l,t) = 0, m + k + 1 = 0 (3.1) 


and 


G 2 { 1 , m, k, t) + G 2 {m, k, I, t) + G 2 {k,l,m,t) = 0 , 1 + m + k = 0 


(3.2) 


Here, we have removed the singularity of A: = 0 and / = 0, motivated by the supposed continuous 
distributions of Gi and G 2 in the wave number space for all time, especially in the limit of 
t +CX 0 . Of course, the above two hold if we take Gi = 0 and G 2 = 0 as the initial conditions 
for the artihcial transient phase before the asymptotic states emerge. 

Next, from a numerical simulation perspective and the consideration of simplicity, there is 
the necessity for the adoption of a hnite support for Di, ZI 2 , and D 4 , and consequently hnite 
supports for G 2 and Ukk following from fl2.57p and fl2.59p : Let maxA; denote the adopted constant 
upper bound for k, which makes it convenient and less time-consuming the computation of the 
integrals involved in those constraints formed in the physical space and the objective. To 
normalize the hnite supports with maxA; and to remove u from the asymptotic state solutions, 
we introduce the non-dimensional scalings through 


t = 


V (maxA:)^ ’ 


k = maxA;k, Ukk{k,t) = 




(k,i) 


max/c 


G 2 { 1 , m, k, t) = 


G^\l, rh, k, i) 
(max/c)^ 


Di{m, k, I, t) 


D3{m,k,l,t) 


1^4 £){a) 

(max/c)^ 

1^4 k,i,i) 

(max/c)"^ 


D2{m,k,l,t) 


A;, /, t) 


1^4 

(max/c)^ 

1^4 

(maxA;)® 


(3.3) 


which is an extension of fl2.90p . It then follows that the correlations in the physical space are 
scaled through 


r =--, Uij(r,t) = iu maxA:)^ A), WiuAr^t) = z/^(maxA;)^tCja;/“Vf, A), 

maxA; 

bJiUJj{Y, A) = i/^(maxA;)'^cuiaJj'^“^(f, A), Uijk{r,s,t) = {v maxA;)^ A), 

U{ij)KL{r, s. A) = (z/ maxA;)'^ s. A), Q{y, A) = (z/ maxA;)"^ 

Qi^L(r, s,A) = (z/ maxA;)'‘Q^^(f, s,A), Wi{y.)wj{yi)uk{y)oJi{'z) = z/'‘(maxA;)® WiWjOJkOJi^°^ 


(f,s,A), 
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a;i(x)wj(x)wfc(y)^^Kz) = y^{Yi\ayikf ■ ■ ■ 

Under the above scalings, fl2.57p and fl2.59p are, respectively, transformed into 


(3.4) 


y 

Ifi 


+ 2k-‘W£{k,t) 


pi /•min(l,/+fc) 

= 47r / di d|k + i|/|k + i| (1-02) 

Jo j\i-k\ 


X 


2+li!h!±iihU!UGr>(ik+ii./.M') 


0 


|k + l|2 

fc/(l-02) 
Ik + ip 


iG^^\\k+\\,k,l,i) 


|k+l|2 -/2 

2 ki 


, 0 = 


(3.5) 


and 


^+p+/ 2 ^ ik+ir 

dt 


|k + i|,fc,t) 


U)$“^(|k + i|,fc,/,t) -D^“^(|k + i|,/,fc,t) (3.6) 


All the constraints listed in Subsection 12.51 and constraint fl3.2p maintain the same structures 
after the transformation - the mere replacement of the original dimensional quantities with the 
corresponding dimensionless. The dehnition of support itself implies that 

D^°'\k,m,i,i) = G2°'^(fc, m,/, t) = 0, fc > 1 or m > 1 or / > 1; Ukk{k,i) = 0, k>l (3.7) 


Here, the inclusion of > 1 is to account for the operations like k —k — 1, rri = k + l^ —k, 
etc. The objective function in fl2.87p takes the equivalent form of 



dkk^Lrll}(k,i) 


to be maximized at each and every great i 


(3.8) 


We have thus a transformed SOCP problem for the asymptotic states with the magnitude of 
its wave-numbers bounded by 1 and parameter-free. 

Thirdly, under the above scaling, fl2.88p becomes 


k(7l“)(0,() = -Stt p dkk*ui^(kj) < 0 (3.9) 

The question is whether the linear forms of fl3.5p and fl3.6p and the decaying behavior of fl3.9p 
result in the automatic satisfaction of the variance inequalities fl2.79p at large time in the 
asymptotic states. 

Computationally, we need to integrate fl3.5p and fl3.6p starting from i = 0 with appropriately 
prescribed initial conditions for and that satishes the constraints fl2.6ip and fl2.30p . 
and there is an artihcial transient phase before the emergence of the asymptotic states. 
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3.1 Sub-model 


For the numerical simulation of the asymptotic states, we need to discretize the unit cube 
support = [0,1]^ for -0^“^ and the constraints and fl3.5p and fl3.6p . 

and we obtain a large-scale SOCP problem. To have a rough estimate about the number of 
discrete variables involved, we take a uniform mesh size of 1/50 along each edge of the cube 
which results in 51^ nodes for the mesh, and accordingly, there are 51^ x 5 discrete variables 
associated with and at each discretized time instant. Meanwhile, there are a large 

number of the discretized constraints to be imposed. Based on this estimate and the computing 
resources available, we will focus on the asymptotic states of the sub-model which involves only 
the second and the third order correlations, fl2.30p . fl2.6ip . fl2.68p through fl2.70p . fl3.2p . 


and (I3.8p . which are recorded below. 


Ji 


+ 2P|C/“(i,i) 


nl /•min(l,/-|-Ai) 

■ 4:71 di —0^) 

Jo J\t-k\ 


X 


2pQ^ + kie-P\ f ^ . 

2 H-- ) fc G 2 (m, I, k, t) 


+ 0 + 


fcf(l-02) 




k, I, i) 


, 0 = 


— k'^ — r 

2 ki 


(3.10) 


> 0 


(3.11) 


<4 


U<f (s' -T, i) - u!f (s-r, t) - ujf (s', i) + ujf (s, i) 

+ 2 a (Ujf (0, i) - U^f (s' - s, ()) - 2 ,3 (0, i) - Ujf (r, i)) 

- a 13 (u^f (s' - f, i) -Utf(s-T, i) - ulf (s', i) + ulf (s, i)) 

df (0, i) - df (f. + «' (df (0> - df (§' - 1)) 

+«(di*(“' ~ “ df (s - f, i) - df *)+df (s, i) 

X yf (o,i)- df (S' - s, i)+/?ddf (0. f - df (f. f 


/3(df(S'-i-.l)-df(s-f.«)-df(s'.*) + df(s.<)) . y = 11.12 (3.12) 


r(“) / 




(«)/ 


(s' — f, f) — tCja;/“^(s — f, f) — WiU)/°'\s', i) + i) 
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+ apywiUjj 


,(a)(g' _Y,i)— WiUj^°'\s — r, t) — Wiu:j^°'\s\ i) + Wiu/'^\s, i)j 


<4 


X 


ujf (0, i) - ulf (f, i) + a'^ (plf (0, i) - uP (s' - s, 

+ a {uP (s' - f , () - up (s - r, () - (s', i) + (7/f (s, 

(0, i) — uJjUJj^^'^ (s' — s, t) + /3^ fuJjUJj^^^ (0, i) — ojjoj/^^ (r, 


— f3 (s' — r,t} — ujjOi}/‘^\s — r, t) — ujjUj^°‘\s', i) + t) j 


, tj = 11,12 

(3.13) 


(s' — f, t) — (s — f, t) — (s', t) + (s, i) 


+ 2 Q;^ci;jCi;/“^(0, t) — — s, t) j —2/9 ^ci;ja;/“^(0, t) — 0 Ji 0 Jj^°‘\r, t) j 

upjjrio) (^s _ f ^ (s', t) + (s. 


<4 


(0, t) - (r, i) + (0, t) - (s' - s, t) j 


X 


+ a (s' - f, i) - cuicu/") (s - f, i) - cuicu/") (s', i) + (s, 

oJjuFj^^^ (0, i) — 0 JjUj^°^ (s' — s, t) + (^Uju/^'^ (0, t) — (r, t) j 

(“) (^g' — Y,i) — (s — r,i) — (s', i) + (s. 


UjUj 


tj = 11,12 

(3.14) 


(j' 2 °'^(m, k, I, i) + G'' 2 \ 1 , k,m,i) = 0, rh + k + 1 = 0 


Pa) 


(3.15) 


{l,rh,k,i) + G) 2 \'rh,k,l,i) + G) 2 \k,l,m,i) = ^, 1 + rh + k = 0 (3.16) 


pa) , 


pa) 


G 2 ‘^^(fc,'m,/, t) = 0, fc>lorm>lor/>l 


(3.17) 


and 


G^ 2 \k, rh, i, i) 


< 1 


(3.18) 


dk P f7^fc^(/c, t) to be maximized at each and every great i 


(3.19) 
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Here, the addition of fl3.18p is to preserve the scaling property of fl2.89l) for the sub-model, 
due to the non-enforceability of fl2.79p : The upper bound is normalized to the unity, allowed 
by the scaling property of the equations involved; This bound restricts the feasible domain, 
and thus, it plays the role of reducing computing time. The above is a SOCP problem when 
discretized. We should mention that within the sub-model, the equality constraint fl2.27p is 
satisfied automatically which lends the basis to take (13.1 Op . Also, the general property, fl2.88p 
or (13.9p . holds within the sub-model. 

The constraints of (I3.15p and (13.161) can be simplihed for the sake of numerical simulation. 
Under (13.151) . the equality constraint (I3.16P reduces to and is equivalent to 

0 ^ 2 "^ (/, m, k, i) + 0 ^ 2 ^ (m, k, f, i) -|- (fc, /, m, t) = 0, k <i < m <k + i, m = —k — i 

(3.20) 

To prove this claim, we first infer from (13.151) the following two specihc representations com¬ 
patible with (I3.17p . 

fc,/, f)-|-G 2 “^(/, fc, m, t) = 0, |fc —/| < m < min|fcl}, 0<fc,/<l (3.21) 

and 

G'" 2 \i,k,i,i) = 0, 0 < fc < min|2/, 1 }, 0 < / < 1; 

G 2 °'^(m, fc,/, t)-|-G 2 ‘^^(/, fc, m, t) = 0, m — / < < min|ml}, 0</<m<l (3.22) 

using \k — l\ < |k -|- 1| < k + I and the like. Next, we verify the equivalence between (I3.16P and 
(I3.20p under (13.151) . (I3.2ip and (I3.22p in the three cases oirh = I, m > I and m < I, respectively. 
The case, rh = f, is trivial in that (I3.16p is satished by (13.2ip and (I3.22p i ii k > 1 and by (I3.22p 
if A: < /. The case, m > I, includes three subcases: k = I, k < I and k > 1. For k = I, (I3.16p 
is satished by (I3.22p . For k < I, (I3.16P is satished by (I3.20p . Subcase k > I contains two 
possibilities of m > k > I and k > rh > 1: the hrst possibility is equivalent to Subcase k < I, 
which can be demonstrated with (I3.15P and (I3.22p : the second possibility is equivalent to the 
hrst. The last case, m < /, is equivalent to the case of rh > /, which can be verihed with the 
application of (I3.15P to (I3.16p . 

We can argue for the equivalence between (I3.2ip and (I3.22p based on that (I3.22p can be 
rewritten in the equivalent form of 

G 2 “^(rh, fc,/, f)-|-G 2 “^(/, fc, rh, t) = 0, rh — / < < min|rhl}, 0</, rh<l (3.23) 

The equivalence between (I3.15P and (I3.2ip comes from that the latter is a specihc representation 
of the former and the latter implies the former in that for any rh, k and 1 with rh -|- k -(-1 = 0 
involved in (I3.15p . ih = —(k -f 1) and rh = |k -l- 1| G [\k — l\, k + 1], there exist a corresponding 
set of {k, l,rh} in (I3.2ip . Consequently, (I3.15P and (I3.22p are equivalent. 

The above results imply that we can replace (I3.15P and (I3.16P with (I3.20p and (I3.22p in 
numerical simulations. 

At every hxed time instant t, each of the constraints, (I3.12p through (I3.14p . involves a scalar 
function of the components of r, s, s' and their diherences and the non-negative constants a 
and f3 to be hxed. To make the SOCP problem computationally feasible, under specihcally 
hxed a and /9, an adequate set of the collocation points in the physical space has to be selected 
at which the constraints are to be imposed. 
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3.2 Discretization Scheme 

To carry out the numerical simulation, we discretize the unit cubic support [0,1]^ with a 
structured hexahedral mesh of uniform size 6k, 


JV —i JV —1 1 \ —1 

'G^-)= U U IJ 


N-l N -1 N -1 

T>^(a) = 

2 - - - 

m=i n 2 =l ^ 3=1 

77,2,ns) = [kn^,kn^ + 5k) x [kn^,kr, 

The half-open intervals are adopted to avoid double couuLUig, m uuc cvaiuanuii ui <^2 i 
interval of each axis should be closed to represent fully the closed support, which is not crucial 
in practice due to the enforcement of (13.171) . 

The distribution of k, I, i) in cubic cell 'H{ni, n 2 , ns) is approximated by 

distribution. 


^712 + ^^) ^ [^ns) kn^ 6k) (3.24) 

are adopted to avoid double counting in the evaluation of 0 ^ 2 '^; the end 

Tiilrl ViP plnQprl tn rpr^rpcipnt fnllv tVip plnciprl cmr^nnrt ^xrViir 


■ the tri-linear 


G' 


/c,/, f;ni,n2,ns) 


_^( 0 ) / . I . I T +\ ^ ^ ^ 

— Go (ni -h 1, n2 -f 1, ns + 1, t)-^-^-^— 

6k 6k 6k 

-y(O), 


I ^(0),^ I 0 I 1 ^ ^*^1 ^«'2 + l ^ ^ ^713 

G 2 (ni -|- 1, n2, ns -1-1, t) -^-^-^— 

6k 6k 6k 

I /^{O)/■ ,3 I 1 ^n-i+l — TTl k — kn 2 1 ~ kn^ 

G 2 (t^i, n2 4- 1, ns -l- 1, t) -^-^-^— 

6k 6k 6k 

. ^(0) / I H 2 \ ^rii+l kh kn 2 +l k i kn 3 

+ G^2 [ni, n2, ns + 1 , t) ---^ 

6k 6k 6k 

I ,^(0) ( ,3 ,3 ~ k — kn2 kn^+l — 1 

-\- G 2 (ni -h 1 , n2 -h 1 , ns, t) -^-^-^- 

6k 6k 6k 

^(0) / ... ^ ^^2+1 ^ ^ns+l I 

+ G 2 + ^,n2,n3,t) -^-^-^- 

6k 6k 6k 

I /^(O)/ I 1 ^ni+l ifl k kn2 kn3+l ^ 

+ G^ (ni, n 2 + 1 , ns, t) ---- 

6k 6k 6k 

+ Gf (m, n2, ns, i) kn3+i-l 

6k 6k 6k 

where G 2 °^(ni,n 2 ,ns,f) denotes the value of G^*^^ at node (ni,n 2 ,n 3 ) at time instant i. The 
distribution of G^\rh, k, 1, i) in V (a) is approximated by 

A^-1 

G2“^(m, k,i,i) = E X[k )(^) XlK2,K2+i)ik) X[fc„3, fc„3+i)(0 <^2 k, i, i; m, n2, ns) 

711,122,713 = 1 

(3.26) 
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where X[kn^,kn^+i) the like are the characteristic functions. 

The direct constraints for ''^ 2 ,'n.a, t) can be established by the substitution of fl3.25p 

into fl3.17p . fl3.18p . fl3.20p and fl3.22p . First, the boundary conditions of the hnite support (13.17p 
require that 

G 2 °^(ni, 77 - 2 , n-s, t) = 0, rii = N or n 2 = N or = N (3.27) 

Next, we evaluate the specihc representation of symmetry (13.221) at the nodes and the mid¬ 
points between the nodes to obtain 

G® (rii, 772, n-i, t) = 0, 1 < 772 < min{2 77i, iV — 1}, 1 < ^71 < iV — 1; 

G^°^(77i,772,773,t) G^°^ ( 773 , 772, 77i, t) = 0, 

^73 — 77i < 772 < min{77i -|- 773 , — 1}, 1 < 77i < 773 < — 1 (3.28) 

Thirdly, the evaluations of the asymptotic state constraint (13.201) at the nodes and the mid¬ 
points between the nodes result in 

G2°^ (771,772, 773 , f) -h G2°^(772,773,77i,t) G^^^ ( 773 , 77i, 772, t) = 0, 

1 < 771 < 772 -I- 1, 772 < 773 + 1, 77,3 < min{77i -|- 772, A^ — 1} (3.29) 


Finally, the bound constraint (I3.18p leads to 


G® (771, 772 , 773, t) 


< 1 


(3.30) 


To implement constraints (13.111) through (I3.14p . we need to represent in terms of G 2 
by solving (13.101) . To this end, we integrate (I3.10p in the time interval re [t, t -|- (5t] to get 

ut>(k,i+si) 

— exp ( - 2 P St) Ulf (k, i) 

pi pmin{lJ-\-k) 


+ 4:71 dl 

Jo j\i-k\ 


dm lrh (l — 0^) 


X 


- 2 + 


2/202 + ^f0_/2 


777^ 


] ^ J [2 P (r — t — ] 0^2'^ (m, /, k, t) 


0 + —-^rrTT —if dr exp [2 P (r — t — 6i) ] G^^^ {m, kj^r) 


777 ^ 


0 = 


777 ,^ — P — 
2ki 


(3.31) 


we then approximate G^\m, k, I, r) in a linear fashion. 


G 2 °^ (tti, k, i, t ) = (777, fc, /, i) —- -f G 2 °^ (777, k,i,i + 6i) '^ ^ , re \i,i + 6t\ (3.32) 

St St 


29 





































and integrate dr analytically to obtain 


u'£(k,i + si) 

— exp{—2k^6t)Ulf{k,t) 

TT 1 — {1 + 2 k'^ 6i) exp(—2 k'^ 6i) 


16 6t 


dl dm ( ^ . f. -2 f2 r2 - 2 ?2 - 2 

X / — / - (m^ + k^ + - 2 P f - 2 P - 2 f rri^ 

Jo I J\i-k\ m V 

X (TTi'^ — k'^ — + 2 k"^ P') 0 ^ 2 '^ (m, k, I, i) 


m* + P + 3 P — 2 P 2 ^ 2 “^ (m, 1, k, i) 


TT 1 — 2 P (5t — exp(—2 P 5t) 


16 




P d/ /•min(l,[+fc) ^2 - 2 ?2 - 2 

X / — / — f + k +l -2k I -2k m -2rm 


'o I J\i-~k\ 


X 


m 


— P — P + 2 P P ) gP (m, k, I, i + (5t) 


— (+ P — P P + 3 P — 2 P 2 (m, /, p t + (5t) 


(3.33) 


Next, substituting fl3.26p into fl3.33p and integrating dm analytically (if desired) 

I L — rC I j/Ctj, ^ j 

result in 


Ut^{k,i + 5t) 

— exp{-2k‘^Si)uji‘l}{k,t) 

TT 1 — (1 + 2 P (ft) exp(—2 P (ft) 

16 (ft P 

iV-l 1 

X E E G^°^(ni + mi, n 2 + m 2 , Ug + mg, t) Mmpm 2 m 3 {ni,n 2 , Ug; p 

ni,n2,n3=l mi,m2,m3=0 

TT 1 — 2 P 5t — exp(—2 P 5t) 

^ 1^ P 

iV-l 1 

X E E G 2 °^(ni + mi,n 2 + m 2 ,ng + mg,t + (fp M,nim.2mz{ni,n2,no] k) 

ni,n2,77-3 = 1 7711,7712,7713=0 

(3.34) 
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Here, each Mmirn 2 m 3 {ni,n 2 ,n^]k) involves corresponding 1-dimensional integrals dl and 

kn2 

jfc„ 3 +i 2 -dimensional integrals (if the above-mentioned analytical integration not imple- 

kng 

mented) 


kn2+l ^ ^min{l,Z+fc,fcni+l} 

dl / dTTL, 

Jma.x{\i—k\,knj^} 


ng+i ^ ^min{l,/+fc,fc„3+i} 

dl / dm 

5 J vaax{\i-k\,kn-^} 


The details of Mmim 2 m 3 {ni,n 2 ,n 3 ] k) are not given here dne to the cnmbersomeness and they 
can be easily inferred from the procednre listed above. 

We can now select a hnite set of collocation points for k G (0,1) to enforce fl3.1ip . 

U^^kikm^i + 5i) m = (3.35) 

For the evalnation of fl3.12p throngh fl3.14p . we resort to the straight-forwardly derived 

^Ulf{y,i + 5i) 


= / dk Fij{T, k) exp{—2 k'^ di) {k,i) 


IT 


N-l 


16 6i 


E E Gf\ni mi, n2 + m 2 , -F m 3 , i) 


ni,n2,n3=l mi,m2,m3=0 


11 1 ^ r\ 1 “ (1 + 2 di) exp(—2 k'^di) 

X / dkFijyV^kj ]^mim2m^ (?7,i,n2,n3; k) 




TT 


Af-1 


E E G®(ni mi, 77,2 + m 2,773 + m3,t -1- di) 


16 di 


ni,n 2 ,n, 3 =l mi,m2,7713=0 


X / dkFij{r,k) 

Jo 


1 — 2k‘^ di — exp (—2 di) 

¥ 


M 


mim2m3 


>1,772,^3; k) 


(3.36) 


where 


, /sin(ffc) sin(ffc) cos{rk)\ (f^P — 3) sin(ffc)-|-3ffccos(ffc) 

rij[T,k) :=0jj/c| -;--h 


/^3^2 


(3.37) 


With the help of (I3.34p . the objective fnnction of (I3.19p can be recast, at f -|- di, in the 
eqnivalent form of 


N-l 1 

E E G 2 °^( 77 i -h mi, 772 + rn 2 ,n 3 -f- m 3 ,i + di) 

ni,n2,n3=l mi,m2,m3=0 


¥ jT 1 ~ 2 P hf — exp (—2 P ht) 

X / dk -'^mim2m3 (^l5 ^2? ^3; (3.38) 

Jo k^ 
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To carry out the numerical simulation of the discretized SOCP problem above, we need to 
start from t = 0 and fix the time step 6t. We notice the inconsistency between this t = 0 and 
the supposed asymptotic states at large time; Computationally, we adopt appropriate initial 
conditions for 0) and G 2 °^(ui, U 2 , na, 0) satisfying the discretized constraints above, solve 

for G® (ui, 77-2, Us, 6i) and 6t} through the optimization, sequentially in time, and expect 

the asymptotic state solutions to emerge from these artihcial transient states. 

The dimensionless turbulent energy ai t + 6t is 


u‘£(o,i+si) 


= 4:71 I dk exp{—2 Si) {k,i) 


71 

4 si 


N-l 


Gf\ni + mi, 712 + m 2 ,712, + m2„i) 


ni,n 2 ,n 3 =l m\,m2,mz=0 


1 "^ 1 — (1 + 2 P (5t) exp(—2 (5t) - 

^ / dk r- -^mim2m3 ('^l; 712, 712, k) 


TT 


N-l 


k^ 


4 St 


E E Gf\7ii + mi, 712 + m2,7i3 + m 3 ,i + Si) 


ni,722,^3=1 mi,m2,m3=0 

^ - 1 — 2 k"^ si — exp(—2 Si) 


X dk 


k^ 


■Sdmxm2mzi+li k) (3.39) 


which can be computed with the known ul^{k,i) and G^\ni,7i2,n2,i) and the newly solved 

G2°^(?7,i,?7,2,n3,t + si). 


References 

[1] CVX Research, Inc., CVX: Matlab software for disciplined convex programming, version 
2.0. http://cvxr.com/cvx, April 2011. 

[2] P. A. Davidson, Turbuleiice qti iTitroductioTi for Scientists and Engineers, Oxford University 
Press, New York, 2004. 

[3] M. C. Grant and S. P. Boyd, The CVX Users’ Guide Release 2.1. CVX Research, Inc. 
October 24, 2014. 

[4] T. Hahn, CUBA-a library for multidimensional numerical integration. Computer Physics 
Communications 168 (2005) 78-95. 

[5] T. Hahn, The CUBA library. Nuclear Instruments & Methods in Physics Research A559 
(2006) 273-277. 

[6] http://www.feynarts.de/cuba 


32 







[7] M. S. Lobo, L. Vandenderghe, S. Boyd and H. Lebret, Applications of second-order cone 
programming. Linear Algebra and its Applications 284 (1998) 193-228. 

[8] E. E. O’Brien and G. C. Francis, A consequence of the zero fourth cumulant approximation, 
J. Fluid Mech. 13 (1962) 369-382. 

[9] Y. Ogura, A consequence of the zero fourth cumulant approximation in the decay of 
isotropic turbulence, J. Fluid Mech. 16 (1963) 33-40. 

[10] I. Proudman and W. H. Reid, On the decay of a normally distributed and homogeneous 
turbulent velocity field. Phil. Trans. Roy. Soc. London. Series A 247 (1954) 163-189. 

[11] L. Tao and M. Raniakrishna, Multi-scale turbulence modelling and maximum information 
principle. Part 1. arXiv:1009.1691vl [physics.flu-dyn] 2010 . 

[12] L. Tao and M. Raniakrishna, Multi-scale turbulence modelling and maximum information 
principle. Part 2. arXiv:1307.4888v3 [physics.flu-dyn] 7 Jan 2014 . 

[13] L. Tao, Multi-scale turbulence modelling and maximum information principle. Part 3. 
arXiv:1408.0376v2 [physics.flu-dyn] 15 Mar 2015 . 

[14] T. Tatsumi, The theory of decay process of incompressible isotropic turbulence. Proc. Roy. 
Soc. London. Series A 239 (1957) 16-45. 


33 





